#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include "noise.h"
#include "random.h"

#define TABSIZE 0x100
#define TABMASK 0xff

#define CURVE(t) (t * t * (3.0 - 2.0 * t))
#define LERP(t, a, b) (a + t * (b - a))
#define INDEX(x) (perm[(x) & TABMASK])
#define INDEX2(x, y) (perm[((x) + perm[(y) & TABMASK]) & TABMASK])
#define INDEX3(x, y, z) (perm[((x) + perm[((y) + perm[(z) & TABMASK]) & TABMASK]) & TABMASK])

static int perm[TABSIZE];
static float vt1[TABSIZE];
static float vt2[TABSIZE][2];
static float vt3[TABSIZE][3];

float pcg_vnoise1(float x) {
	float t, v0, v1;
	int xi;
	
	xi = floor(x);
	t = x - xi;
	
	t = CURVE(t);
	v0 = vt1[INDEX(xi)];
	v1 = vt1[INDEX(xi + 1)];
	
	return LERP(t, v0, v1);
}

float pcg_vnoise2(float x, float y) {
	float t, v00, v01, v10, v11, r1, r2;
	int xi, yi;
	
	xi = floor(x);
	yi = floor(y);
	
	v00 = vt1[INDEX2(xi, yi)];
	v01 = vt1[INDEX2(xi + 1, yi)];
	v10 = vt1[INDEX2(xi, yi + 1)];
	v11 = vt1[INDEX2(xi + 1, yi + 1)];
	
	t = x - xi;
	t = CURVE(t);
	r1 = LERP(t, v00, v01);
	r2 = LERP(t, v10, v11);
	
	t = y - yi;
	t = CURVE(t);
	return LERP(t, r1, r2);
}

float pcg_vnoise3(float x, float y, float z) {
	float tx, ty, tz, v00, v01, v10, v11, r1, r2, r3;
	int xi, yi, zi;
	
	xi = floor(x);
	yi = floor(y);
	zi = floor(z);
	
	v00 = vt1[INDEX3(xi, yi, zi)];
	v01 = vt1[INDEX3(xi + 1, yi, zi)];
	v10 = vt1[INDEX3(xi, yi + 1, zi)];
	v11 = vt1[INDEX3(xi + 1, yi + 1, zi)];
	
	tx = x - xi;
	tx = CURVE(tx);
	r1 = LERP(tx, v00, v01);
	r2 = LERP(tx, v10, v11);
	
	ty = y - yi;
	ty = CURVE(ty);
	r3 = LERP(ty, r1, r2);
	
	v00 = vt1[INDEX3(xi, yi, zi + 1)];
	v01 = vt1[INDEX3(xi + 1, yi, zi + 1)];
	v10 = vt1[INDEX3(xi, yi + 1, zi + 1)];
	v11 = vt1[INDEX3(xi + 1, yi + 1, zi + 1)];
	
	r1 = LERP(tx, v00, v01);
	r2 = LERP(tx, v10, v11);
	
	r1 = LERP(ty, r1, r2);
	
	tz = z - zi;
	tz = CURVE(tz);
	return LERP(tz, r3, r1);
}

float pcg_gnoise1(float x) {
	float t, v0, v1, dx, dx1;
	int xi;
	
	xi = floor(x);
	
	dx = x - xi;
	dx1 = dx - 1.0;
	
	v0 = vt1[INDEX(xi)] * dx;
	v1 = vt1[INDEX(xi + 1)] * dx1;
	
	t = CURVE(dx);
	return LERP(t, v0, v1);
}

float pcg_gnoise2(float x, float y) {
	float t, v00, v01, v10, v11, r1, r2, dx, dy, dx1, dy1;
	int xi, yi, i;
	
	xi = floor(x);
	yi = floor(y);
	
	dx = x - xi;
	dy = y - yi;
	dx1 = dx - 1.0;
	dy1 = dy - 1.0;
	
	i = INDEX2(xi, yi);
	v00 = vt2[i][0] * dx + vt2[i][1] * dy;
	i = INDEX2(xi + 1, yi);
	v01 = vt2[i][0] * dx1 + vt2[i][1] * dy;
	i = INDEX2(xi, yi + 1);
	v10 = vt2[i][0] * dx + vt2[i][1] * dy1;
	i = INDEX2(xi + 1, yi + 1);
	v11 = vt2[i][0] * dx1 + vt2[i][1] * dy1;
	
	t = CURVE(dx);
	r1 = LERP(t, v00, v01);
	r2 = LERP(t, v10, v11);
	
	t = CURVE(dy);
	return LERP(t, r1, r2);
}

float pcg_gnoise3(float x, float y, float z) {
	float tx, ty, tz, v00, v01, v10, v11, r1, r2, r3, dx, dy, dz, dx1, dy1, dz1;
	int xi, yi, zi, i;
	
	xi = floor(x);
	yi = floor(y);
	zi = floor(z);
	
	dx = x - xi;
	dy = y - yi;
	dz = z - zi;
	dx1 = dx - 1.0;
	dy1 = dy - 1.0;
	dz1 = dz - 1.0;
	
	i = INDEX3(xi, yi, zi);
	v00 = vt3[i][0] * dx + vt3[i][1] * dy + vt3[i][2] * dz;
	i = INDEX3(xi + 1, yi, zi);
	v01 = vt3[i][0] * dx1 + vt3[i][1] * dy + vt3[i][2] * dz;
	i = INDEX3(xi, yi + 1, zi);
	v10 = vt3[i][0] * dx + vt3[i][1] * dy1 + vt3[i][2] * dz;
	i = INDEX3(xi + 1, yi + 1, zi);
	v11 = vt3[i][0] * dx1 + vt3[i][1] * dy1 + vt3[i][2] * dz;
	
	tx = CURVE(dx);
	r1 = LERP(tx, v00, v01);
	r2 = LERP(tx, v10, v11);
	
	ty = CURVE(dy);
	r3 = LERP(ty, r1, r2);

	i = INDEX3(xi, yi, zi + 1);
	v00 = vt3[i][0] * dx + vt3[i][1] * dy + vt3[i][2] * dz1;
	i = INDEX3(xi + 1, yi, zi + 1);
	v01 = vt3[i][0] * dx1 + vt3[i][1] * dy + vt3[i][2] * dz1;
	i = INDEX3(xi, yi + 1, zi + 1);
	v10 = vt3[i][0] * dx + vt3[i][1] * dy1 + vt3[i][2] * dz1;
	i = INDEX3(xi + 1, yi + 1, zi + 1);
	v11 = vt3[i][0] * dx1 + vt3[i][1] * dy1 + vt3[i][2] * dz1;
	
	r1 = LERP(tx, v00, v01);
	r2 = LERP(tx, v10, v11);
	
	r1 = LERP(ty, r1, r2);
	
	tz = CURVE(dz);
	return LERP(tz, r3, r1);
}

float pcg_fBm1(float x, pcg_noise1_fn_t prim, float H, float lacunarity, float octaves) {
	float value, rem;
	int i, oi;
	
	value = 0.0;
	oi = (int) octaves;
	for (i = 0; i < oi; i++) {
		value += prim(x) * pow(lacunarity, -H*i);
		x *= lacunarity;
	}
	
	rem = octaves - oi;
	if (rem > 0.0) 
		value += rem * prim(x) * pow(lacunarity, -H*i);
	
	return value;
}


float pcg_fBm2(float x, float y, pcg_noise2_fn_t prim, float H, float lacunarity, float octaves) {
	float value, rem;
	int i, oi;
	
	value = 0.0;
	oi = (int) octaves;
	for (i = 0; i < oi; i++) {
		value += prim(x, y) * pow(lacunarity, -H*i);
		x *= lacunarity;
		y *= lacunarity;
	}
	
	rem = octaves - oi;
	if (rem > 0.0) 
		value += rem * prim(x, y) * pow(lacunarity, -H*i);
	
	return value;
}

float pcg_fBm3(float x, float y, float z, pcg_noise3_fn_t prim, float H, float lacunarity, float octaves) {
	float value, rem;
	int i, oi;
	
	value = 0.0;
	oi = (int) octaves;
	for (i = 0; i < oi; i++) {
		value += prim(x, y, z) * pow(lacunarity, -H*i);
		x *= lacunarity;
		y *= lacunarity;
		z *= lacunarity;
	}
	
	rem = octaves - oi;
	if (rem > 0.0) 
		value += rem * prim(x, y, z) * pow(lacunarity, -H*i);
	
	return value;
}

float pcg_multifractal1(float x, pcg_noise1_fn_t prim, float H, float lacunarity, int octaves, float offset) {
	float value;
	int i;
	
	value = 1.0;
	for (i = 0; i < octaves; i++) {
		value *= (prim(x) + offset) * pow(lacunarity, -H*i);
		x *= lacunarity;
	}
	
	return value;
}

float pcg_multifractal2(float x, float y, pcg_noise2_fn_t prim, float H, float lacunarity, int octaves, float offset) {
	float value;
	int i;
	
	value = 1.0;
	for (i = 0; i < octaves; i++) {
		value *= (prim(x, y) + offset) * pow(lacunarity, -H*i);
		x *= lacunarity;
		y *= lacunarity;
	}
	
	return value;
}

float pcg_multifractal3(float x, float y, float z, pcg_noise3_fn_t prim, float H, float lacunarity, int octaves, float offset) {
	float value;
	int i;
	
	value = 1.0;
	for (i = 0; i < octaves; i++) {
		value *= (prim(x, y, z) + offset) * pow(lacunarity, -H*i);
		x *= lacunarity;
		y *= lacunarity;
		z *= lacunarity;
	}
	
	return value;
}

float pcg_hetero1(float x, pcg_noise1_fn_t prim, float H, float lacunarity, float octaves, float offset) {
	int i;
	float value, inc, rem;
	
	value = offset + prim(x);
	x *= lacunarity;
	
	for (i = 1; i < octaves; i++) {
		inc = (prim(x) + offset) * pow(lacunarity, -H*i);
		inc *= value;
		value += inc;
		x *= lacunarity;
	}
	
	rem = octaves - (i - 1);
	if (rem > 0) {
		inc = (prim(x) + offset) * pow(lacunarity, -H*i);
		inc *= value;
		value += inc * rem;
	}
	
	return value;
}

float pcg_hetero2(float x, float y, pcg_noise2_fn_t prim, float H, float lacunarity, float octaves, float offset) {
	int i;
	float value, inc, rem;
	
	value = offset + prim(x,y);
	x *= lacunarity;
	y *= lacunarity;
	
	for (i = 1; i < octaves; i++) {
		inc = (prim(x,y) + offset) * pow(lacunarity, -H*i);
		inc *= value;
		value += inc;
		x *= lacunarity;
		y *= lacunarity;
	}
	
	rem = octaves - (i - 1);
	if (rem > 0) {
		inc = (prim(x,y) + offset) * pow(lacunarity, -H*i);
		inc *= value;
		value += inc * rem;
	}
	
	return value;
}

float pcg_hetero3(float x, float y, float z, pcg_noise3_fn_t prim, float H, float lacunarity, float octaves, float offset) {
	int i;
	float value, inc, rem;
	
	value = offset + prim(x,y,z);
	x *= lacunarity;
	y *= lacunarity;
	z *= lacunarity;
	
	for (i = 1; i < octaves; i++) {
		inc = (prim(x,y,z) + offset) * pow(lacunarity, -H*i);
		inc *= value;
		value += inc;
		x *= lacunarity;
		y *= lacunarity;
		z *= lacunarity;
	}
	
	rem = octaves - (i - 1);
	if (rem > 0) {
		inc = (prim(x,y,z) + offset) * pow(lacunarity, -H*i);
		inc *= value;
		value += inc * rem;
	}
	
	return value;
}

float pcg_distnoise1(float x, pcg_noise1_fn_t prim, pcg_noise1_fn_t dist, float scale) {
	float dx;
	dx = dist(x);
	return prim(x + scale * dx);
}

float pcg_distnoise2(float x, float y, pcg_noise2_fn_t prim, pcg_noise2_fn_t dist, float scale) {
	float dx, dy;
	dx = dist(x,y);
	dy = dist(y,x);
	return prim(x + scale * dx, y + scale * dy);
}

float pcg_distnoise3(float x, float y, float z, pcg_noise3_fn_t prim, pcg_noise3_fn_t dist, float scale) {
	float dx, dy, dz;
	dx = dist(x,y,z);
	dy = dist(z,x,y);
	dz = dist(y,z,x);
	return prim(x + scale * dx, y + scale * dy, z + scale * dz);
}

void pcg_snoise(int s) {
	int i, j, k;
	float rho, theta;
	
	pcg_sdrand(s);
	
	// initialize index array
	for (i = 0; i < TABSIZE; i++) {
		perm[i] = i;
	}
	
	// shuffle index array
	for (i = 0; i < TABSIZE; i++) {
		j = pcg_drand() * TABSIZE;
		k = perm[j];
		perm[j] = perm[i];
		perm[i] = k;
	}
	
	// initalize values
	for (i = 0; i < TABSIZE; i++) {
		vt1[i] = 1.0 - 2.0 * pcg_drand();
		
		theta = pcg_drand() * M_PI * 2.0;
		vt2[i][0] = cosf(theta);
		vt2[i][1] = sinf(theta);
		
		rho = pcg_drand() * M_PI;
		theta = pcg_drand() * M_PI * 2.0;
		vt3[i][0] = cosf(theta) * sinf(rho);
		vt3[i][1] = sinf(theta) * sinf(rho);
		vt3[i][2] = cosf(rho);
	}

}